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We present a representative set of analytic stationary state 
solutions of the Nonlinear Schrodinger equation for a sym- 
metric double square well potential for both attractive and 
repulsive nonlinearity. In addition to the usual symmetry 
preserving even and odd states, nonlinearity introduces quite 
exotic symmetry breaking solutions - among them are trains 
of solitons with different number and sizes of density lumps in 
the two wells. We use the symmetry breaking localized solu- 
tions to form macroscopic quantum superposition states and 
explore a simple model for the exponentially small tunneling 
splitting. 



I. INTRODUCTION 

Many features of Bose-Einstein condensates (BECs) of 
dilute atomic gases in a single well external potential at 
zero temperature are well described by mean field the- 
ory [1,2]. In the mean field picture all condensate atoms 
have the same macroscopic wave function satisfying the 
Gross-Pitaevskii (GP) equation. In this paper we investi- 
gate the stationary states of BEC in a symmetric double 
square well potential. We find analytic solutions of the 
GP equation for both symmetry preserving and symme- 
try breaking stationary states of the attractive and repul- 
sive nonlinearity. The solutions presented in the paper 
give such analytic expressions for what are seen to be 
stationary soliton trains in the double well - among them 
are such trains with different number and sizes of density 
lumps in the two wells. Single dark solitons [3,4], bright 
soliton [5] and soliton trains [6] have been recently experi- 
mentally observed in trapped BECs, suggesting that their 
double well analogs may be experimentally accessible. In 
addition we present, as an application of the mean field 
symmetry breaking solutions, a zero order macroscopic 
mean field descriptions of macroscopic quantum super- 
position states (Schrodinger Cat state) in a double well 
BEC system. 

Symmetry breaking mean field solutions, such as we 
observe in this exact treatments, are expected in the at- 
tractive case as an attractive condensate in the ground 
state tends to localize in one well or the other. Symmetry 
breaking solutions for a nonlinear Schrodinger equation 
was first pointed out in the context of molecular states [7] . 
Symmetry breaking mean field states for repulsive con- 
densates have been discussed in the two-state model of 



condensate dynamics in a double well [8-11], and seen in 
the nonlinear numerical studies of the GP equation in a 
symmetric quartic double well [12]. The present analytic 
work thus confirms the numerical work of D'Agosta and 
Presilla in Ref. [12] in the context of a double square well. 
Such macroscopic quantum self-trapped states have also 
appeared on the studies of transport on a dimer modeled 
by discrete nonlinear Schrodinger equation [13]. 

BECs in a double well and multi-well systems have 
been studied in the context of coherence [14], Josepson 
tunneling [8,15,16], squeezed states [17], the superfluid to 
Mott transition [18] and condensate fragmentation [19]. 
In discussions of condensate tunneling it is well known 
that a high barrier leads to condensate fragmentation 
in which two or more distinct single particle states are 
macroscopically occupied. For a repulsive condensate, 
raising the barrier leads to the condensate in the two 
wells from being coherent to being incoherent in a Fock 
state [19]. The analysis herein gives the nonlinear modes 
of the entire double well in a mean field picture when 
all the atoms have the same single particle wavefunction. 
Correlation effects leading to condensate fragmentation 
are neglected here and thus the theory presented applies 
directly only to the case of strong tunneling. However, 
the mean field states obtained could form the basis for a 
correlated description. 

The GP equation is a cubic nonlinear Schrodinger 
equation(NLSE) [20] where the particle interactions give 
rise to such effective nonlinearity. The NLSE has 
been successful in modeling many other natural phe- 
nomenon besides BEC. It describes light pulses in op- 
tical fibers [21], helical excitations of a vortex line [22], 
Bose-condensed photons [23] , spin waves in magnetic ma- 
terials [24], and disordered media [25]. Despite being 
a canonical physics problem [26], the symmetric dou- 
ble square well problem has not, to our knowledge, been 
solved for nonlinear Schrodinger equation. Although the 
discussions in the paper is exclusively for Bose-Einstein 
condensates, the analysis will apply to any system satis- 
fying cubic NLSE. 

The symmetry breaking localized one particle mean 
field states can be used to form a zero order two- 
configuration Schrodinger cat states of the form 4>%,f t ± 
bright- There have been several reports of the creation 
of Schrodinger cat states in various condensed matter 
systems [27,28]. In the context of BEC, several au- 
thors have suggested producing such states [11,29-32], 
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although none have been demonstrated experimentally. 
In a double well, as is found analytically in this paper, 
the mean field ground state for an attactive condensate is 
a symmetry breaking state localized in one of the wells. 
The superposition of such degenerate localized states is 
a "cat" state. We calculate the tunneling splittings for 
such states using correct mean field single particle states 
starting from the full N-body Hamiltonian. Such two- 
configuration tunneling splittings are exponentially small 
in the N-body wave function overlap. Particle correla- 
tions are still neglected, but strong mean field effects ac- 
counted for. 

The article is organized as follows. In Sec. II we present 
the full set of symmetry preserving and symmetry break- 
ing analytic solutions of stationary NLSE for a symmetric 
double square well potential. In Sec. Ill we discuss an ap- 
plication of the symmetry breaking solutions - the possi- 
bility of creating superpositions of macroscopic quantum 
states, and calculate the tunneling splittings of such cat 
states. Remarks and discussions in Sec. IV conclude the 
paper. 

II. DOUBLE SQUARE WELL 

The stationary NLSE with a potential has the form 

[-dl + V | f{x) | 2 +V tra "{x) } f{x) = n f{x) , (1) 

where fix) is the mean field condensate wavefunction 
in the longitudinal direction, /x is the eigenvalue or the 
chemical potential, and ij is the nonlinearity parameter 
which is proportional to the number of atoms and the 
s-wave scattering length. All quantities in Eq.(l) are di- 
mensionless. 

Analytic solutions of the GP equation for harmonic 
and quartic double well potentials are not possible, so we 
have chosen to investigate the infinite square well with 
symmetrically placed finite rectangular potential barrier. 
The potential is of the form 

!oo, \x\ > a 
0, b<\x\<a (2) 
V , \x\ < b 

For clarity, Fig.l shows a picture of this potential. Dou- 
ble well traps can be created in experiments with a com- 
bination of optical and magnetic trapping. Varying the 
laser strengths the barrier or the depth and the width 
of the trap can be easily tailored to experimental speci- 
fications. The double well traps created in experiments 
usually have gaussian barriers; however, the qualitative 
behavior of the stationary states of such wells would be 
the same as discussed in this paper for a double square 
well. 

We present the analytic solutions of Eq.(l) with the 
potential Eq.(2). Solutions in an infinite well and a finite 
well have been presented for both attractive and repul- 
sive condensates [33-35]. In Eq.(l) r\ > corresponds 



to repulsive condensate while 77 < corresponds to at- 
tractive condensate. The solutions of NLSE in a zero 
potential are Jacobian Elliptic functions [36] . Such func- 
tions are well known in the soliton literature, and also 
as the solution to the anharmonic classical oscillator, i.e. 
9 + 9 — 8 3 /3l = 0. An example of the standard nota- 
tion for Jacobi Elliptic functions is sn(a; | to), where to 
is the elliptic parameter. The period is given by 4K(m), 
where Kim) is the complete elliptic integral. The value 
of to is bounded between and 1. It interpolates the 
elliptic functions between trigonometric and hyperbolic 
functions. There are 12 elliptic functions all of which are 
solutions to the NLSE. Of the 12 elliptic functions, six 
are bounded and six are unbounded. Of the six bounded 
functions, only sn(x | m), cn(x | to), dn(x | to) have dis- 
tinct physical forms. Others differ only by a translational 
shift or a rescaling of the amplitude. The six unbounded 
functions can be represented as a quotient of the above 
three functions in different combinations. We will find 
that the pieces of these unbounded functions arc those 
appropriate in the barrier region of the double well for a 
repulsive condensate. Table I summarizes the functions 
relevant to this work. 

Solutions in the three regions will be written in the 
form 

{fiix), -a<x<-b 
Mx), \x\ < b (3) 
foix), b<x <a 

The solutions vanish on and outside the hard wall bound- 
ary at I a; I > a. The solutions will be found subject to 
continuity of fix) and fix) at x = ±6 and the normal- 
ization condition f_ a dx\fix)\ 2 = 1. The vanishing of the 
solutions at the hard walls is taken as built into the ellip- 
tic functions and does not form an additional boundary 
condition. The solutions are divided into two different 
categories - Symmetry preserving and Symmetry break- 
ing. Taking advantage of the symmetry of the problem 
finding symmetry preserving states reduce to solving a 
set of three nonliear algebraic equations. The symmetry 
breaking states require solving five simultaneous nonlin- 
ear equations which is a far more difficult undertaking. 

A. Symmetry preserving states 

Symmetry preserving states are the states that pre- 
serve the symmetry of the N-particle many-body Hamil- 
tonian. Simply put, they are the even and odd solutions. 
As we will find out in the next section, there can also be 
solutions which does not preserve even/odd symmetry 
expected from linear quantum mechanics. 
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1. Attractive nonlinearity 

Symmetric solutions take the following form 

fi(x) = A cn(k(x + a) - K(m) \ m) , (4a) 
f 2 (x) = A 2 dn(k 2 x + K(m 2 ) | m 2 ) , (4b) 
f 3 (x) = A cn(k(x - a) + K(m) \ m) , (4c) 

and antisymmetric solutions take the form 

fi(x) = Acn(k(x + a) — K(m) | m) , (5a) 
f 2 (x) = A 2 cn(k 2 x + K(m 2 ) | m 2 ) , (5b) 
/ 3 (ar) = -Acn(k(x-a)+K(m)\m), (5c) 

where A, A 2 , k, k 2 , m and m 2 are free parameters. 
fi(x) and fs(x) have been chosen to preserve odd and 
even parity. Note that the elliptic parameter K{m 2 ) dis- 
places the an in the barrier region to make it antisymmet- 
ric. In the next section we describe uniquely nonlinear 
type solutions which does not preserve such parity. The 
condition that the states vanish at the hard walls at a 
and -a are built into the form of the solutions. 

Symmetric and antisymmetric solutions are solved us- 
ing the same method. Substituting the symmetric solu- 
tions into Eq. (1) with the potential Eq. (2), following 
conditions are obtained 

A 2 = 2mk 2 /r], A\ = 2k\jr] (6a) 
fi = (1 - 2m)k 2 , n={m 2 - 2)k 2 2 + V a (6b) 

The boundary condition fi(—b) = f 2 (—b) is equivalent 
to f 2 (b) = h(b), and requires 

Acn(kuj - K(m) \ m) = A 2 dn(-fc 2 6 + K(m 2 ) | m 2 ) (7) 

where w = a — b is the width of each of the wells. Conti- 
nuity of the first derivative requires 

Aksn(kuj — K(m) \ m) dn(kuj — K{m) \ m) = A 2 m 2 k 2 
x sn(-fc 2 fe + K(m 2 ) \ m 2 ) cn(-fc 2 6 + K(m 2 ) \ m 2 ) (8) 

Finally, the normalization condition is 

2A 2 ( dxcn 2 (k(x - a) + K(m) \ m) 

Jb 

f -b 

+2A\ dx dn 2 (fc 2 x + K(m 2 ) \ m 2 ) = 1 . (9) 
Jo 

Eq. (9) can be written as 

2 -§ [E(k 2 b + K(m 2 ) | ma) - E(m 2 )} - ^(1 - m)w 

+ 2 -£ [E{m) - E(-ku; + K(m) | m)] = 1 . (10) 

where E(k 2 l \ m) is standard notation for an incomplete 
elliptic integral [36]. 

Equating of Eqs. (6b) gives us a constraint on the en- 
ergy. Substitution of Eqs. (6a) into Eqs. (7), (8), and 



(10) produces a system of four simultaneous equations - 
an energy condition, a nontrivial normalization and two 
enforcing the continuity of the wavefunction and its first 
derivative at the interior discontinuity of the potential. 
The four equations can be reduced to three equations in 
three unknowns. These are 

y/mkcn(ku; — K(m) \ m) 
= Adn(-A6 + K(m 2 ) \ m 2 ) (11a) 

^pmk 2 sw(kuj — K(m) \ m) dn(kuj — K(m) \ m) = m 2 X 2 
sn(-Afr + K(m 2 ) | m 2 ) cn(-Afr + K (m 2 ) | m 2 ) (lib) 

f [E(Xb + K(m 2 ) | ma) - E(m 2 )} - ^ (1 - m)u 

+f[E(m)-E(-kuj + K(m)\m)} = l. (11c) 

where A = \J V "~ 2 -m^ k = ^2 and ui = a — b. This is a 
system of three nonlinear algebraic equations with three 
unknown variables m, m 2 and k and four experimental 
parameters - the box width 2a, barrier height V a , barrier 
width 2b and nonlinearity parameter r\. 

This system of equation Eqs. (11) is analogous to the 
set of equations for linear Schrodinger equation for a par- 
ticle on a box double well potential [26]. However, the 
normalization equation Eq. (11c) here is nontrivial and 
gives an additional condition. These equations can ide- 
ally be solved by a multidimensional secant method, and 
that is the method we use to find the roots. However, the 
nonlinear parameter space is too large to choose a good 
starting point for the roots to converge. As we will see in 
the next section when we deal with a set of five equations 
for the symmetry breaking solutions, it is almost impos- 
sible to find the roots and the analytic solutions without 
a good initial choice of parameters from an approximate 
numerical solution. 

Such numerical approximations to the exact solutions 
of Eq.(l) with the double well potential Eq.(2) can be 
generated by the shooting method [37]. However, the 
cubic nonlinearity generated from the mean-field inter- 
actions of the atoms introduces numerical stiffness into 
the resulting two-point boundary value problem. To ac- 
curately compute the numerical solutions, Gear's meth- 
ods [38] are employed which are efficient in overcoming 
the numerical stiffness by utilizing backward differencing 
formulas. The resulting shooting scheme is then easily 
implemented and both the normalized symmetry preserv- 
ing and symmetry-breaking states are computed along 
with their chemical potential. We note that by adjust- 
ing the shooting angle, the normalization to unity can be 
satisfied. 

Knowing the chemical potential and the value of the 
solution at barrier boundary x — b from the shooting rou- 
tine numerics we can find the three approximate roots of 
Eqs. (11). With the form of the solutions and the ap- 
proximate roots at hand, secant method is used to solve 
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the Eqs. (11) to find the exact analytic solutions. In 
Fig. 2 we show the first four odd and even states. The 
states are ordered according to the chemical potential [i. 
A barrier height of V = 100, barrier width of 2b = 1/5, 
well width 2a = 1 and nonlinearity of r\ = —100 were 
used. Table II shows the solution parameters for Fig. 2. 
The true mean field ground state for an attractive con- 
densate in this case is a symmetry breaking state where 
the condensate localizes in one well or the other as is de- 
scribed in the next section. The first excited even state 
for this well in Fig. 2(c) where the condensate has one of 
the peaks on top of the barrier is a uniquely nonlinear 
state [12] which does not have any counterpart in linear 
Schrodinger equation. For fi > V all even solutions are 
of this kind, however even for fi < V strong nonlinearity 
can give rise to such states. Symmetric solutions of this 
kind has the form f 2 (x) = A 2 cn(k 2 x | to 2 ). 

The antisymmetric solutions were found using a 
similar method. For reference the system of equa- 
tions is Y / m2 _ Acn(— Xb + K(m 2 ),m 2 ) — \/mkcii(k(a — 
b) — K(m), to); y /m 2 X 2 sn(— Xb + K (m 2 ), m 2 )dn(— Xb + 
K(m 2 ) 7 m 2 ) = y/rnk 2 sn(ku; — K(m), m)dn(kuj — 
K(m), to); f [E(Xb + K{m 2 ) | m 2 ) - E(m 2 )} - (1 - 
m)uj-^- (l-m 2 )b+f(E(m)-E(-kw + K(m) | to)) = 

1, where A = ^^'^ = k 2 . 

We would like to note that unlike linear quantum me- 
chanics, for attractive condensate the eigenvalue or chem- 
ical potential of the antisymmetric state for this well di- 
mensions has a lower value than the symmetric case. This 
behaviour is only true for strong nonlinearity The total 
energy per particle for the antisymmetric state is how- 
ever always greater than the symmetric case. Similar 
behavior of symmetric and antisymmetric state chemical 
potentials has also been found in the case of ring poten- 
tials [33]. 

2. Repulsive nonlinearity 

Symmetric solutions take the form 

fi{x) = Asn(k(x + a) | to) , (12a) 
f 2 (x) = A 2 ds(k 2 x + K(m 2 ) | m 2 ), (12b) 
f 3 (x) = -Asn(k(x - a) | to) , (12c) 

and antisymmetric solutions take the form 

fi(x) = Asn(k(x + a) \ m) , (13a) 
f 2 (x) = A 2 cs(k 2 x + K(m 2 ) \ m 2 ) , (13b) 
fz{x) = Afin(k(x — a) | m) , (13c) 

Substitution of theses solutions into Eq. (1) with the dou- 
ble well potential Eq. (2) gives the following equations for 
the amplitude and the chemical potential 

A 2 = 2m,k 2 /r], A\ = 2k\ /r, (14a) 
H = (l + m)fc 2 , [i = -{2m 2 -l)kl + V (14b) 



Just like for the attractive case the three simultaneous 
equations obtained from the boundary conditions, nor- 
malization and the energy conditions are following 

Y / rofcsn(fco; | m) = Ads(— Xb + K{m 2 ) \ m 2 ) , (15a) 

^fmk 2 cn(fcw | m) dn(fcw | m) = —A 2 x 
cs(-A6 + K(m 2 ) | to 2 ) ns(-A6 + K(m 2 ) | m 2 )(15b) 

4X 2 b/r] - 4X 2 m 2 b/r] + ^-uj 
+ f[cs(-Xb + K(m 2 ) | m 2 )dn(-Xb + K(m 2 ) \ m 2 ) 
— cs(Xb + K{m 2 ) | TO 2 )dn(A6 + K{m 2 ) \ m 2 )\ 
-f[-E(Xb + K(m 2 ) | to 2 ) + E(-Xb + K(m 2 ) \ m 2 )] 

-fE(koj | m) = 1 (15c) 

where A = y ^ 1Jr \ l } 2 ^ ~ = k 2 . A similar set of equations 
is obtained for the antisymmetric case. 

The ground state and the first three symmetry pre- 
serving excited states are shown in Fig. (3). The well 
dimensions used here are different than the attractive 
case which was chosen to show the peculiarities of attrac- 
tive condensate. A barrier height of V = 1000, barrier 
width of 2b = 1/10, well width 2a = 1 and nonlinearity 
of rj = 100 were used here. Table II shows the solution 
parameters for Fig. 3. In addition to the even and odd ex- 
cited states there can also be symmetry breaking states 
as described in the next section. For a repulsive conden- 
sate the lowest symmetry preserving state is always the 
ground state. 

B. Symmetry breaking states 

Symmetry breaking states are uniquely nonlinear 
states where different size or number of "lumps" are 
present in the two wells. Such stationary states with 
strong localization and different number of nodes in the 
two symmetric wells are not possible for linear Sturm- 
Liouville systems. Finding such solutions confirms and 
extends the numerical work [12] and the two state tun- 
neling models [8-10,39] of the double well where macro- 
scopic quantum self-trapping has been predicted. On the 
N-particle level the stationary states should preserve the 
symmetry of the Hamiltonian and can only be symmet- 
ric and antisymmetric. So these asymmetric states arise 
due to the nonlinearity associated with the mean field 
approximation. 

In the work of D Agosta and Presilla [12] a non-linear 
trial function and relaxation method for patial differ- 
ential equations was used to numerically find both the 
symmetry preserving and symmetry breaking states of 
the GP equation in a symmetric harmonic/quartic dou- 
ble well. The difficult task of choosing the right trial 
functions and the possibility of false minima leading to 
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artifacts in such methods motivated us to treat the model 
double square well potential and to find the roots of these 
algebraic equations, and thus find the exact analytic so- 
lutions. The qualitative behaviour of solutions in any 
symmetric double potential should be the same as ours, 
and wherever the set of parameters we used overlaps with 
those of Ref [12] there is one-to-one correspondence in the 
solutions. 

1. Attractive nonlinearity 

Solutions with no nodes inside the barrier region take 
the form 

fi(x) = A\ cn(fci(x + a) — K(m\) \ mi) , (16a) 
h(x) = A 2 dn(k 2 (x + d) + K(m 2 ) | m 2 ) , (16b) 
f 3 (x) = A 3 cn(k 3 (x - a) + K(m 3 ) | m 3 ) , (16c) 

and solutions with nodes inside the barrier are 

fi(x) = A\ cn(fci(a; + a) — K(m\) \ mi) , (17a) 
f 2 (x) = A 2 cn(k 2 {x + d) + K(m 2 ) \ m 2 ) , (17b) 
f 3 (x) = -A 3 cn(k 3 (x - a) + K(m 3 ) | m 3 ) , (17c) 

d in Eqs. (16) and (17) is a measure of how far the solu- 
tion under the barrier is displaced from being symmetric. 
The amplitudes and the chemical potentials are 

A\ = 2mikl/rf, A\ = 2k 2 2 jm A\ = 2m 3 fc§/»7 (18a) 
H = (l-2mi)fe?, 11= (m 2 -2)k 2 + V , 
pi = (l-2m 3 )fc| (18b) 

The set of five equations in five unknowns are 

y^mT a cn(A 3 (mi) | mi) = /3dn(Ai(d, m 2 ) \ m 2 ), (19) 

% A^3 7cn(A4(m 3 ) | m 3 ) = (3 dn(A 2 (<i, m 2 ) \ m 2 ) , (20) 

y/m[ a 2 sn(A 3 (mi) | mi) dn(A 3 (mi) | mi) 
= mlP 2 sn(\i(d,m 2 ) | m 2 )cn(Ai(d, m 2 ) \ 2 ) , (21) 

v / m 3 "7 2 sn(A 4 (m 3 ) | m 3 )dn(A 4 (m 3 ) | m 3 ) 
= m^/3 2 sn(A 2 (d, m 2 ) | m 2 ) cn(A 2 (d, m 2 ) | 2 ) , (22) 

- ^1(1 - m 3 )cj + %[E(m 3 ) - E(\,(m 3 ) \ m 3 )\ 

- 2si(l - mi)w + f[E(mi) + E(X 3 ( mi ) \ mi)} 
f[E(X 2 (d,m 2 ) | ma) - B(Ai (d,m 2 ) | m 2 )] = 1 (23) 

where « = ^T^r = k u p = v /gJ = fc 2 , 

T = \J i-2m 3 = fc 3, Ai(d,m 2 ) = fc 2 (rf - 6) + if(m 2 ), 
A 2 (d, m 2 ) = fc 2 (d + 6) + K(m 2 ), A 3 (mi) = aw - if (mi) 



and A 4 (m 3 ) = —71^ + K{m 3 ). This is a set of five non- 
linear equations in five unknowns mi, m 2 , m 3 , d and /x. 
A similar set of equations are obtained for the solutions 
that has nodes inside the barrier. 

As described in the previous section we use a shoot- 
ing method to find the approximate numerical solutions. 
Knowing the eigenvalue and the values of the functions 
at the barrier boundaries at x = ±6, we can reduce five 
equations with five unknowns to equations with two un- 
knowns. With just two unknowns we can use a graphical 
method [35] to find the approximate solutions. Such ap- 
proximate roots are then used to find the exact analytic 
roots of these five equations using a multidimensional se- 
cant method, and thus we obtain the analytic solution of 
the symmetry breaking states. Without a good bracket- 
ing on the roots obtained from first solving it numerically 
its extremely unlikely for a secant method to converge for 
a set of five nonlinear equations. 

For Fig. 4 we use a well dimension of 2a = 1, 2b = 1/10, 
V = 1000 and nonlinearity r\ = —100. We use a differ- 
ent well dimension here than the attractive symmetric 
case just to show the varieties of asymmetric solutions. 
Table III shows the solution parameters for Fig. 4. The 
solutions can be classified as multiple node solutions - 
zero node, one node, two node and such. The lowest 
symmetry breaking state for attractive condensate is the 
ground state of the system as the clustering of parti- 
cles in one of the wells minimizes the energy for strong 
enough self interaction. There can be ground and ex- 
cited state solutions with assymetrically placed peaks on 
top of the barrier. The analytic form of such solutions 
is f 2 (x) = A 2 cn(k 2 (x + d) I m 2 ). For an asymmetric 
ground state, increasing the barrier height localizes the 
condensate more into the well, on the other hand increas- 
ing the barrier width pushes the peak of the condensate 
density more towards the center of the well on top of the 
barrier. 

Fig. 2 and Fig. 4 are the bright soliton solutions in a 
double well. It shows the one, two, three and four soliton 
solutions. Bright soliton and soliton trains have recently 
been observed in attractive condensates of 7 Li [5,6]. Un- 
like stationary soliton trains of equal density lumps in a 
single potential well, double well geometry has stationary 
soliton train solutions with unequal density lumps as is 
shown in Fig. 4. There exists a whole class of such many- 
soliton solutions. As an example, Fig. 5 shows an ana- 
lytic solution of a symmetry breaking cight-soliton bright 
soliton train in a well of dimensions 2a = 1, 2b = 0.1, 
V = 1000, and for nonlinearity 77 = —500. 

2. Repulsive nonlinearity 

Solutions with no nodes inside the barrier are 

fi{x) = Aisnihix + a) \m x ), (24a) 

f 2 (x) = A 2 ds(k 2 (x + d) + K(m 2 )\m 2 ), (24b) 

f 3 (x) = A 3 sn(k 3 (x - a) | m 3 ) , (24c) 
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Solutions with nodes inside the barrier are 

fi(x) = Ai sn(fci(x + a) \ mi) , (25a) 
f 2 (x) = A 2 cs(k 2 (x + d)+K(m 2 ) | m 2 ), (25b) 
f 3 (x) = -A 3 sn(k 3 (x - a) | m 3 ) , (25c) 

The amplitude and chemical potentials for the states that 
has no nodes inside the barrier are 

A\ = 2m 2 k\jr\,A 2 2 = 2kl/rj,A\ = 2m 3 k\^ (26a) 

Hi = (1 + mi)kl,fi 2 = (1 - 2m 2 )k\ + V a , 

H 3 = (l+m 3 )k 2 3 (26b) 

For reference the equations are 

Y / m7asn(a uo \ mi) = /3ds(Ai(d, m 2 ) \ m 2 ) (27) 

v / mi'7sn(-7w | m 3 ) = f3ds(\ 2 (d,m 2 ) \ m 2 ) (28) 

v / rrira 2 cn(aw, mi)dn(aw, mi) 
= -/3 2 ns(Ai(d, m 2 ) | m 2 )cs(\i(d, m 2 ) \ m 2 ) (29) 

A /m 3 7 2 cn(— -fiv | m 3 )dn(— -fu \ m 3 ) 
= -/3 2 ns(A 2 (d, to 2 ) \ m 2 )cs(X 2 (d,m 2 ) \ m 2 ) (30) 

4/3 2 6/r/ - 4f3 2 m 2 b/r] + ^-uo + ^f-u 

+ ^[cs(Ai(d,m 2 ) | TO 2 )dn(Ai(rf,TO 2 ) | m 2 ) 
- cs(\ 2 (d, m 2 ) | m 2 )dn(A 2 (d, m 2 ) \ m 2 )] 
- fE(au; | mi) - f E(ju; \ m 3 ) 

- f[-E(X 1 (d, ma) | ma) + E(X 2 (d, m 2 ) | m 2 )] = 1 (31) 

where the same notations as in the attractive case has 
been used. Here the five unknown variables are mi, m 2 , 
m 3 , d, and fi; a, (3 and 7 are functions of the elliptic 
parameters and the chemical potential /i. A similar set 
of equations is obtained for the solutions that has nodes 
inside the barrier. 

The first four states are plotted in Fig. 6 for a non- 
linearity of r\ = 100 and the same well dimension as the 
repulsive symmetry preserving case, 2a = 1, 2b = 1/10 
and V = 1000. Table III shows the solution parameters 
for Fig. 6. Again the solutions can be classified as one- 
node, two-node, three-node symmetry breaking states. 
For repulsive condensates the asymmetric ground state 
has a much higher energy and is in fact the second ex- 
cited state of the double well. Note that for the two 
two-node solutions keeping one node inside the barrier 
and another outise the barrier is energetically more fa- 
vorable than having two nodes outside the barrier. In 
Fig. 7 we show the symmetry breaking ground state as 
we change the nonlinearity. It evolves from being almost 
localized for small nonlinearity to having three distinct 
density lumps for high enough nonlinearity. 



III. SCHRODINGER CAT STATE OF BEC IN A 
DOUBLE WELL 

As was shown in the previous section, the mean field 
ground state of attractive condensate and some of the 
excited states of both attractive and repulsive conden- 
sate are symmetry breaking states. For the symmetry 
breaking localized states such as the attractive ground 
state, coherent quantum tunneling between the degener- 
ate states removes the degeneracy and forms a superpo- 
sition of the mean field states. Such localized superpo- 
sition states of the form 4>f e j t ± bright are Schrodingcr 
cat states. On the other hand, the usual even and 
odd symmetry preserving delocalized states of the form 
^ N = (<pieft ± 4*right) N are not traditional Schrodingcr 
cat states. For the cat states, tunneling splitting is ex- 
ponentially small in the N-body wave function overlap. 
In the following we find the zero order two-configuration 
mean field cat state tunneling splitting starting with the 
N-particle Hamiltonian with pseudopotential interaction. 
It has not gone unnoticed that the ground state of the 
attractive condensate is cat-like [40,41]. Cirac et al. [11] 
have studied the ground state of Josephson-coupled two- 
species condensates which has similarities with conden- 
sate in a double well. In the next section we deliberate 
on the experimental realization of cat states of BEC in a 
double well. 



A. Schrodinger cat state tunneling splitting 

The N-body Hamiltonian for a system of N weakly 
interacting identical bosons each of mass m in an external 
potential V ext is 

N / h 2 \ 

H » = E ( + Wro) + 1/2 E^ r ^-)' 

(32) 

Here T^r^r.,) = gS(ri — r^) is the Fermi 'contact' 

pseudo-potential, and g = 4nc ^ h where a s is the s-wave 
scattering length characterizing the binary atomic colli- 
sions. 

For a fully condensed Bose condensate the N-body 
wavefunction can be written as a symmetric product of 
single-particle wave functions 

<Mri,r 2 , . . . , r N ) = 0(ri)0(r 2 ) . . . <j>(r N ) = <f> N (33) 

where </>(i"i)'s are the single particle mean field wavefunc- 
tions normalized to unity J dr\(f>\ 2 = 1. 

The expectation value gives us the N particle mean 
field energy 

(<P N \H N \<p N ) =»N- N{N 2 +1) g J dr\^ (34) 
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where /i is the chemical potential. We can generalize 
these to the left and right localized GP solutions in a 
double well, which in the above equatin would corre- 
spond to replacing <p by <Pl and <f>R. The expectation 
value with respect to the left and right localized states 
contains overlap integrals which no longer vanishes be- 
cause of their non-orthogonality - 



N(N-l) 



(A) 1 



-gN 2 J drfflfafto (Af- 1 + [iN (A) N (35) 



where A = J dr<j)* L (j>ji is the overlap integral. The even 
and odd combinations of the left and right localized so- 
lutions <j>[ift ^ bright are a two configuration model for 
Schrodinger cat superposition states. Taking <f>L and <j)R 
to be real, the expectation value of the energy at this 
simplest level of approximation is 



Es.A — 



H N \<ft)± 



H 



nVPr 



N ) 



l±{<t>L\<t>R) 



N 



(36) 



Although this equation is identical in appearance with 
Eq.(20) of Cirac et al. [11], our use and inclusion of the 
exact mean held effect on the fully localized left and right 
well solutions differ from their treatment of spinor con- 
densates. The tunneling splitting is the difference in an- 
tisymmetric and symmetric energy 



A_E = Ea — Eg 



(37) 



For the case when the overlap is extrememly small and 
for a large number of particles the normalization factor 
in the denominator can be ignored and the splitting can 
be written as 



\N-1 



AE « -2fiN(A) N + 2gN 2 J dr^R^M^Y 

-N(N-l)g f dv^HAf- 2 (38) 



This shows explicitly how the cat state tunneling split- 
ting depends on the overlap of the localized single parti- 
cle mean field wave functions. However, in our calcula- 
tions we find the exact splitting by use of Eq. (37) since 
physically realizable splitting can only be generated for 
a significant overlap such that we cannot completely ig- 
nore the the A N term in the denominator. Since A is 
always less than 1 the splitting is exponentially small in 
the wavefunction overlap. 

We use the solution of the one-dimensional GP equa- 
tion Eq. (1) to find the tunneling splitting and its 
dependence on other quantities. The conversion fac- 
tor to get the energy from a dimensionless quantity is 
fi 2 /(2ml 2 ) [35], where I is the length of the box. To find 
the splittings in one dimension, the coupling constant 
'g' should be replaced by the dimensionless effective one- 
dimensional coupling constant g e ff- The dimensionless 



nonlinearity ij of the NLSE Eq. (1) is related to g e ff by 
the realtionship r\ = g e ffN, where N is the total num- 
ber of particles. For experimental purposes where a con- 
densate is three dimensional or can be quasi-one dimen- 
sional the effective coupling constant g e f / depends on the 
transverse dimensions of the trap, the species of atoms 
(whether attractive or repulsive) and the total number of 
particles in a nonlinear and nontrivial way. Even with- 
out knowing the exact g e f / for realistic three dimensional 
condensate we can explore the dependence of the tunnel- 
ing splittings on the number of particles N and on the ef- 
fective coupling constant. The relationship between the 
effective coupling constant g e // and the transverse di- 
mensions of realistic double well traps that will give the 
correct experimental predictions is under investigation. 



B. Discussions 

Pairs of symmetry breaking mean field states in a 
double well are shown in Fig. 8, coherent tunneling be- 
tween these will produce a cat state. Experimentally 
such macroscopic cat states could be observed by start- 
ing with a localized attractive condensate in the lower 
well of an asymmetric double well potential, and then 
varying the symmetry of the two wells. In Fig. 9 we 
show the log of tunneling splitting for a condensate of 
7 Li as a function of particle number for a double well 
of dimensions 12.5 /i separated by 75 fj, in a box width 
of 100 /j, and barrier height of V = 133. A constant 
effective coupling constant of g e ff = —0.145 has been 
assumed. In Fig. 10 we show the log of tunneling split- 
ting in the same well as a function of \g e ff\ for a fixed 
number of particles - in this case for 500 particles. For 
a cat state, with the addition of more and more parti- 
cles the single particle overlap becomes smaller, and the 
tunneling splitting becomes vanishingly small due to its 
exponential dependence on the overlap and the number of 
particles. Fig. 11 shows the GP single particle tunneling 
splitting between the attractive antisymmetric and sym- 
metric state for g e // = —0.911 which sharply contrasts 
with the cat state tunneling splitting. 

For an example of a cat state, for N = 440 and 
g e ff = —0.145 the peaks of the degenerate states are 
asymmetrically placed on top of the barrier and the sep- 
aration of the peaks is 1.5 /i , the tunneling splitting 48 Hz 
and the tunneling time 21 ms which are within the exper- 
imental range of detection. For higher peak separations 
the overlap is small and the splitting becomes negligblc. 
An optimal cat state with gaussian barriers as is often 
used in experiments where the peaks are well separated 
and the splitting is within the range of detection should 
be attainable with externally tuning the coupling con- 
stant through Fcshbach resonance [42]. The number of 
particles in our study is limited to the order of hundred 
atoms which is within the range of stability of attrac- 
tive condensates [43] such as 7 Li or 85 Rb. Changing the 
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scattering length by Feshbach resonance will allow stable 
attractive condensates to be prepared with several thou- 
sand atoms [6]. For a repulsive condensate, cat states 
may also be prepared making use of the excited localized 
condensate which must be tuned to the right regime to 
get a well localized condensate as shown in Fig. 7(a). 
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IV. CONCLUSION 



We have presented the stationary states of nonlinear 
Schrodinger equation in one dimension for a symmetric 
double square well potential for both attractive and re- 
pulsive nonlinearity. In addition to the symmetry pre- 
serving even and odd states, we find analytic expressions 
for symmetry breaking states that have different num- 
bers and sizes of density lumps in the two wells. For 
attractive condensates these provide the analytic solu- 
tions of the stationary bright soliton trains in a double 
well. Symmetry breaking states do not preserve the even 
and odd parity of the N-particle many-body Hamilto- 
nian. Finding such analytical solution of continous GP 
equation puts the self trapping states as found numeri- 
cally [12], in the 'two-state' tunneling models [8-11], and 
in the discrete nonlinear Schrodinger equation [13] on an 
exact footing. Such unique symmetry breaking states, 
which arc not possible for a linear Scrodinger equation, 
results from the nonlinearity introduced by the mean field 
approximation. 

The superposition of mean field localized states of the 
form (j>g ft ± 

bright are Schrodinger cat states that arise 
due to coherent tunneling between the two degenerate 
states strongly localized in two different wells. Attrac- 
tive condensate in the ground state or repulsive conden- 
sate in its symmetry breaking excited state can be used 
to produce such cat states. In a zero order two con- 
figuration model the splitting is exponentially small in 
the N-body wavefunction overlap. Tailoring the width 
and barrier height of a double well and with adequate 
number of particles in the trap to give the optimal split- 
ting, macroscopic superposition states should be attain- 
able with current BEC technology. 

The use of mean field picture in describing BEC fully 
delocalized in a double well is valid only when the con- 
densate in the two wells are fully coherent. For suffi- 
ciently low tunneling, condensate in a double well cannot 
maintain its coherence and therefore mean field analysis 
of a fully coherent condensate as was presumed here is 
not adequate. Such fragmented condensate with number 
sqeczed configurations can only be treated using theories 
which go beyond mean field theory. However, the avail- 
ability of the mean field analytic solutions as presented in 
this paper provides the zeroth order nonlinear wavefunc- 
tions needed to include important and large mean field 
effects in models which treat fragmentated condensates. 
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FIG. 1. Symmetric double square well potential: the model 
used in this paper 
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FIG. 2. Shown are the first four symmetry preserving 
states for attractive nonlinearity. The barrier walls are at 
x=± 0.1 
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FIG. 3. Shown are the first four symmetry preserving 
states for repulsive nonlinearity. The barrier walls are at x=± 
0.05 
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FIG. 4. Shown are the first four zero-node, one-node and 
two-node symmetry breaking states for attractive nonlinear- 
ity. The barrier walls are at x=± 0.05 
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FIG. 6. Shown are the first four one-node, two-node and 
three-node symmetry breaking states for repulsive nonlinear- 
ity. The barrier walls are at x=± 0.05 
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FIG. 9. Cat State tunneling splitting as described in the 
two configuration model: it shows the exponential dependence 
of the splitting on the number of particles for a fixed coupling 
constant. Energy is in frequency units of Hertz. 
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FIG. 10. Cat State tunneling splitting as described in the 
two configuration model: it shows the exponential dependence 
of the splitting on the effective coupling constant \g e ff\- En- 
ergy is in frequency units of Hertz. 
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FIG. 7. Symmetry breaking repulsive ground state as a 
function of nonlinearity . The barrier walls are at x=± 0.05. 
(a) r7=15, (b) ?7=30, (c) ??=50, (d) r?=100 
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FIG. 11. GP single particle energy splitting between the 
lowest antisymmetric and symmetric state of an attractive 
condensate. The splitting of mean field delocalized states 
slowly increases with particle number, and this runs in a di- 
rection opposite to that of the Cat state tunneling splitting. 
Energy is in frequency units of Hertz 



TABLE III. Solutions parameters for symmetry breaking 
states of attractive and repulsive nonlinearity for Fig. 4 and 
Fig.6. The numbers shown are of sufficient precision as initial 
estimate to be used in the numerical solution of the nonlin- 
ear equations of section II. However, as m — > 1 use of high 
precision arithmetic is required. 





mi 


rri2 


ra 3 


d 


A 4 


FigAa 


0.9999 


1 - 10~ s 


1 - 10" ie 


-0.0680 


-625.27 


FigAb 


1 - 10~ 8 


0.9999 


0.7640 


0.0618 


-174.10 


FigAc 


0.8257 


0.9994 


0.6171 


0.0219 


-62.829 


FigAd 


0.8401 


0.9992 


0.6177 


0.0184 


-62.820 


Fig.&a 


0.9273 


0.9958 


0.2612 


-0.0035 


243.16 


Fig.&b 


0.9273 


0.9959 


0.2529 


-0.0043 


248.95 


Fig.Qc 


0.9612 


0.9992 


0.0016 


-0.0491 


308.14 


Fig.6d 


0.0787 


0.9876 


0.6012 


0.0122 


412.30 



TABLE I. Limits of Jacobian elliptic functions and inte- 
grals. The first two sn and cn are periodic solutions in the 
well while dn, cn, ds, and cs are solutins in the barrier region. 
4K(m) is the periodicity and the elliptic integrals K(m) and 
E(m) both play a role in the system of equations which de- 
scribe the solutions. 



m = m = 1 

sn(tt | m) sin(u) tanh(w) 

cn(u | m) cos(u) sech(tt) 

dn(w | m) 1 sech(tt) 

ds(u | m) csc(u) csch(it) 

cs(w | m) cot(u) csch(u) 

K(m) 7r/2 oo 

E(m) ty/2 1 



TABLE II. Solutions parameters for symmetry preserving 
states of attractive and repulsive nonlinearity for Fig. 2 and 
Fig. 3. The numbers shown are of sufficient precision as initial 
estimate to be used in the numerical solution of the nonlin- 
ear equations of section II. However, as m — > 1 use of high 
precision arithmetic is required. 





m 


m 2 


k 


M 


Fig.2a 


0.9684 


0.9959 


13.25 


-164.42 


Fig. 2b 


0.9758 


0.9935 


13.04 


-161.90 


Fig.2c 


0.6352 


0.9298 


12.47 


-42.03 


Fig.2d 


0.4763 


0.7426 


15.36 


11.18 


Fig.3a 


0.8539 


0.9976 


9.88 


181.06 


Fig.3b 


0.8514 


0.9977 


9.98 


184.51 


Fig.3c 


0.4338 


0.9912 


14.79 


313.75 


Fig.3d 


0.4313 


0.9909 


15.00 


322.24 
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